Influence of boundary conditions on the dynamics of pattern forming systems: the 

case of oscillatory media. 
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■ We study the pattern dynamics in a reaction diffusion model of the activator-inhibitor type 
in the osciUatory regime. We consider finite systems with partially absorptive boundary conditions 

, analizing examples in different geometries in one and two dimensions. We observe that the boundary 

■ conditions have important effects in the pattern forming properties of the systems. In all the 
studied cases, the arising complex behaviour is found te be dependent on the absorption parameter. 
By changing this parameter we can control the asymptotic behaviour to be stationary, periodic, 
quasiperiodic or chaotic. 

o 

CD ■ I. INTRODUCTION 

There are many properties in pattern forming reaction-diffusion systems that do not depend on the geometry of 
. the system (for instance Turing wavelengths, frequencies of oscillation, shapes and velocities of waves or fronts.) 
c/3 101 . Due to this fact, most studies on pattern formation have not paid too much attention to features such as the 
shapes and sizes of the reactors or the boundary conditions, working mainly with infinite systems or periodic boundary 
conditions. However, it is well known that boundaries and boundary conditions can play a relevant role in determining 
the behaviour of this kind of systems. For instance, in reference Q the importance of geometrical characteristics is 
quiet remarkable. Also, recent theoretical papers have shown that the domain shape may have an important effect in 
the propagation of waves in excitable media , and that the consideration of non trivial boundary conditions [Q in 
Q some reaction diffusion-systems leads to qualitatively different dynamic behaviour from the one obtained in infinite 
O [ media. 

Here we analyze problems where the boundary conditions (bc's), as well as the size of the system, are relevant 
parameters determining the kind of complex dynamics arising. We study an activator-inhibitor model in the oscillatory 
regime considering different geometries with partially absorptive bc's. For a density field p(r, t) of a given reactant 
T-H this condition is written as 
^ . 

O. -D{Vp)\^-n^kp\,, (1) 

O 

where D is the diffusion constant, the subscript a indicates evaluation at the boundary, ft is the unit vector normal 
to the boundary in the outer direction, and fc(> 0) is the absorption or albedo parameter. This relation expresses the 
fact that there is an outgoing flux of reactants through the boundary which is proportional to the density. The limit 
k = corresponds to no flux or Neumann bc's, which means a perfectly reflecting (non absorbing) boundary, whereas 
clearly, an increasing value of k corresponds to an increasing absorption at the boundary, and a reduced reflectivity. 
^ , We consider the reaction-diffusion model given by the equations 

; M = V^U + /(m)-U 

O. v = DyV'^v + u--fv-c, (2) 

^ , where the non linear function f{u) = —(u — uq)^ + u characterizes the autocatalytical properties of the activator 

■ u. The model shares the general features of more realistic ones describing chemical systems such as the Belousov 
, Zabotinskii or CIMA reactions For uq = c = 0, it was used in the analysis of thepropagation of fronts in 
' bistable systems and also in the study of pattern formation in inhomogeneous media [^,^. The constant c is here 

set equal to [mo(1 — 7)]. With this election uq represents a translation of the related spatially independent system 
(i.e. Eqs. (|^) without diffusion terms) in the (u, 7;)-planc along the line {u = v). For 7 > 1 the equations describe a 
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bistable medium, while for 7 < 1 they describe an oscillatory one. Here, we consider 7 = .9. In order to work with 
positive values of the densities u and v, we set uq = .4. This last fact is desirable taking into account the kind of bc's 
we will work with, however, it is not strictly necessary, as we have found the same kind of results (in all the problems 
here discussed) using uq — 0. 

For the chosen values of parameters, the homogeneous system associated to Eqs. (|]) has a limit cycle solution 
of period tq ~ 14.66. As initial conditions for our extended system we consider homogeneous states belonging to 
this limit cycle. This (natural) choice is taken in order to focus our attention on the effects caused by the boundary 
conditions. The consideration of inhomogeneous initial states would produce additional pattern forming phenomena 
that would interact with those coming from the boundary conditions. An important characteristic of the extended 
system is that, for £)„ > 1,73, a Turing instability [|l|;0 of wavelength At — 13 coexists (and competes) with the 
Hopf instability ||l|,|I^. 

The aim of this paper is to show that a wide spectrum of possibilities concerning pattern dynamics arises when 
different boundary conditions are considered in an oscillatory medium. We would like to attract the interest of 
experimentalists to study these situations. The experimental observation of the phenomena here predicted might be 
realized by adequately designing chemical systems sharing the properties of our models. Also, the partially absorptive 
conditions may appear in real systems as defects at the boundaries, such as, small separations between the gel and 
the wall of the reactor containing the solution with uncontrolled concentration of reactants. Another experimental 



possibility is to work with electrical systems |11 



All numerical calculations were done as follows. Firstly, the different systems of partial differential equations 
have been approximated by systems of coupled ordinary differential equations, obtained by finite difference schemes. 
Secondly the resulting equations have been solved by a Runge-Kutta 2 method. Different space and time discretization 
schemes were employed in order to check the results. 

The organization of the paper is as follows: in Section II we study a unidimensional system with Neumann boundary 
conditions at one border and partially absorptive at the other. In Section III we study two different bidimensional 
systems corresponding to a circular reactor with partially absorptive boundary conditions and a long rectangular 
system with different boundary conditions on its walls. The last Section is devoted to discussing the results and 
drawing some conclusions. 



II. UNIDIMENSIONAL SYSTEM. 

We now consider Eqs.(^) in the one dimensional domain < a; < L, with Neumann bc's at x = and partially 
absorptive bc's aX x = L: 

dMO) = D.dMO) = 

dxu{L) = -kuu[L) Dydxv{L) -kyv{L). (3) 

As initial conditions we consider homogeneous states belonging to the limit cycle of the associated non extended 
system (Eqs. (^ without diffusion). We analyze the dynamics for different values of L,Dy, fc„ and ky. 

At short times, near x = 0, where there are Neumann bc's, the system evolves very close to the natural homogeneous 
limit cycle of the medium; whereas, near x = L, the motion is perturbed by the albedo bc's which cause the local 
oscillations to be slower and of smaller amplitude. This effect increases with the absorption parameter. Hence, the 
dynamics of the system begins to be inhomogeneous and the perturbation is transmitted from x — L io the rest of 
the system, which asymptotically reaches different dynamical regimes depending on the parameters. The resulting 
regime is independent of the initial condition (when taken on the limit cycle as mentioned above) except maybe in 
some specific parameter regions, as we will discuss later. 

For very high absorption, the oscillations near x = L are inhibited. As a consequence of this, for small enough 
systems (typically L < 2 At, but depending on Dy and fc's), the oscillations are found to be stopped in the whole 
system and stationary patterns are generated. This is because the stationary tendency of the zone near x — L affects 
the whole system. The effect is increased when Dy grows, since the coupling of the rest of the system to the indicated 
zone becomes stronger. (For a higher value of Z?^, a stationary pattern may appear for lower values of the absorption 
parameter or in a larger system.) 

In general, when this stationary regime does not occur, the oscillations are not inhibited. The system asymptotically 
reaches periodic, quasiperiodic, or chaotic inhomogeneous regimes. However, for large enough Dy (typically Dy > 2.3), 
the system evolves toward a Turing pattern, that is also stationary, but of a completely different origin ]lO| than the 
stationary patterns described above. Turing patterns are generated by means of freezing fronts |l^ which travel from 
X = L (the slower oscillating zone) to the rest of the system. As they advance, they inhibit the temporal oscillations 
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and leave behind the characteristic spatially-periodic structures. The process of formation of a Turing pattern is 
shown in Fig.l. 

In Fig. 2 we show the phase diagram indicating the asymptotic behaviour of a system of length L = 70 ~ 5At 
with ku = ky = k as function of Dy and k. The different regions correspond to Turing patterns (TP), inhomogeneous 
periodic oscillations (PO), inhomogeneous quasiperiodic oscillations (QP) and spatiotemporal chaos (CH). 

In the periodic region (PO), the asymptotic behaviours correspond to regimes of periodic traveling waves from 
X — L (albedo be) to a; = (Neumann be). The global period of the motion results always a little bit larger than 
tq. This is because of the early slowdown of the oscillations occurring near x = L, that is later transmitted to the 
whole system. In Fig. 3 we show a spatiotemporal plot of the u-field corresponding to a time window of a typical 
asymptotic regime in the PO region. 

To distinguish between PO, QP, and CH regimes we consider the succession {t^} of times at which u{x — Q,t) 
reaches a local maximum as function of t, i.e. the times when 

d d'^ 

— M(0,i)|t„=0 and ^u(0,t)|t„<0 (4) 

occur simultaneously. Then we define a new succession {p„ = tn — 

The study of the behaviour of {pn} provides a useful and simple way to determine the character of the dynamics 
(periodic, quasiperiodic or chaotic), as was already observed when dealing with inhomogeneous systems In the 
case of periodic motion the value of p„ converges to a constant Poa which coincides with the period of the global 
motion. In the case of quasiperiodic motion the p^s asymptotically show periodic or quasiperiodic behaviour. In the 
chaotic region the p',^s exhibit highly disordered behaviour. Instead of considering the succession we work with 

the values of Pn plotted as a function of tn- Hence we define a function of time p(t) = pn (strictly defined only for 
the discrete values t = tn)- In Fig. 4 we show the typical plots of the p{t) for the PO, QP and CH regimes. In Fig. 
4. a. the early time behaviour can be observed, in which the system at a: = is performing the natural limit cycle 
oscillation of the medium, until the perturbation coming from x = L arrives slowing sown the motion. It is worth here 
mentioning that we have not observed marginal profiles of p{t) where the determination of the kind of dynamics was 
not clear. However, in the transition zones between phases, the transients are long and it is necessary to observe the 
evolution for a very long time. Another remarkable fact is that disordered behaviour such as that appearing in Fig. 
4.C. -here called chaotic- are chaotic in the sense of showing high sensitivity to initial conditions. We have verified 
this fact in a large number of cases by studying the evolution of the distance between solutions with different (but 
close) initial conditions, in the same way as it was done in ||]. 

In the asymptotic regimes of the QP and CH regions, it is not only waves traveling from x = L to x = that 
appear. There are intervals during which waves travel in the opposite sense, and also, from time to time, a center 
emitting waves in both directions may appear in any part of the system (specially in the CH regime) . These effects 
are observed as defects and dislocations in the spatiotemporal plots of the u-fields. In the QP region, the defects 
appear periodically or quasiperiodically in time, while in the CH region they appear in a random way. In Fig. 5 
we show a spatiotemporal plot corresponding to a chaotic regime, together with the plot of p{t) for the same time 
window. There, it can be observed that the arising of defects is related to the occurrence of maxima in p{t), which 
indicates delays in the oscillations at x = 0. 

In a small region of the transition zone between the PO and TP regions, which appears shadowed in the diagram of 
Fig. 2, the behaviour of the system is complicated. In general, a long chaotic transient appears that, later, eventually 
converges to Turing patterns, periodic or quasiperiodic behaviour. In this zone, the asymptotic behaviour is also 
possibly dependent on the initial condition -even when homogeneous and belonging to the indicated limit cycle. We 
have observed also Turing patterns that are generated by freezing fronts coming from a; = (the Neumann be) instead 
of from X = L (the albedo be). One possible explanation for such peculiarities is that, in this zone of the diagram, 
the Turing instability is well established, but the albedo bc's introduce only a very small perturbation (because of the 
small values of fcs) that it is not always strong enough to excite the unstable spatial mode. 

In the phase diagram of Fig. 2 we have indicated with dotted and dashed lines the zones where the transition 
between periodic and nonperiodic regimes respectively occur for a system with (A;„ = 0,ky = k) and (fc„ = k,ky — 0). 
It can be seen that for (fc„ = 0, fc^, = k) the transition occurs for similar -or even lower- values of k than in the case 
ky = ky = k but, for [ky = k,ky — 0), much larger values of k are needed in order to observe nonperiodic regimes. 
This indicates that the system is more sensitive to changes (or imperfections) in the boundary conditions of the u-field 
than of the u-field. 

In the next section we analyze the dynamics of Eqs. (^ with partially absorptive bc's in two different bidimensional 
geometries. 
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III. BIDIMENSIONAL SYSTEMS. 



Circular system. 

Firstly we consider the activator inhibitor model of Eqs.(^) in a circular system of radius R with albedo be for both 
fields 

dru{R) = -kuu{R), D^drv{R) = ~kyv{R). (5) 

where — + is the radial coordinate. We impose homogeneous initial conditions as in the previous section and 
study numerically the radial equations corresponding to (H) : 

1 d du 
^ I d . dv. 

V = Dy- — {r—) + u- jv - c, (6) 
r or or 

with f(u), c, 7 and uo as before. Analyzing the spatiotemporal plots (r,t) of the solutions and the pnS (now defined 
at r = 0), we have observed the same kinds of behaviour (periodic, quasipcriodic chaotic and Turing) as in the 
unidimensional system of the last section. In Fig. 6 we show the phase diagram for two different values of the radius 
of the system R. 

In the first case, for R ~ 70, there is no CH phase because the system is not large enough to allow the required 
disorder on the dynamics of the oscillating phase in the different parts of the system to be produced. However, 
some chaotic isolated points have been observed in the shadowed zone in the transition to the TP region, where the 
properties of the dynamics are similar to those explained for the one-dimensional problem. For R = 80 we have the 
same phases as in the unidimensional problem. Here, the asymptotic regimes correspond to target patterns, that are 
stationary in the TP region, and of waves traveling from r = i? to r = in the PO region, and traveling in both 
directions in the QP and CH regions. 

In the CH region, the radial symmetry is expected to be broken by small angular perturbations in real systems. 

In both diagrams we show the transition lines from periodic to nonperiodic regimes for the cases (ku = k,ky — 0) 
and {ku = 0, fc^, = k). The observed shift of these lines relative respect to that for k^ = ky = k are similar to those in 
the one dimensional problem. The conclusion that the system is more sensitive to boundary conditions on the u-field 
remains valid in bidimensional systems. 



Band shaped system. 

Now we consider Eqs.(^ in the same parameter region as before, in a long rectangular domain {0 < x < L^, —Ly < 
y < Ly) with Ly << with different boundary conditions on its sides. We fix albedo be at a; = with ky — ky — 10 
(which means a strong absorption or very small reflectivity) and Neumann be at a; = L = 80 (fc„ = fc„ = 0). If 
we now fix Neumann be also at the y = iLy sides, we will have (for homogeneous initial conditions) an effective 
one-dimensional system equivalent to the one analyzed in section II. Instead, here we consider another albedo be at 
y = iLy with a different value of k = k^ = ky. In this case the patterns will not have translational symmetry in the 
y-direction. We fix Dy = 2 and study numerically the solutions for different values of k and Ly. The situation to 
analyze is sketched in Fig. 7. 

The one dimensional problem with fc = 10 at one end, and fc = at the other is chaotic for L = 80, Dy = 2. Now, 
in the bidimensional system considered, it is expected that when increasing the albedo parameter on the y = iLy 
sides from zero, the dynamics will begin to differ from the one of the one-dimensional system. In particular, if Ly is 
small and k large, we expect that the system will become stationary because of the stationary tendency introduced 
by this new boundary condition that will constrain the oscillations from the sides. This is in fact what simulations 
have shown. There are also intermediate values of fc and Ly where the stationary tendency of the sides is not enough 
to completely inhibit the oscillations, but it is enough to change their characteristics from chaotic to quasipcriodic or 
periodic. 

In Fig. 8 we show a diagram of the (fc, Lj^)-plane indicating the stationary and non-stationary regions as well as 
some points where periodic, quasipcriodic and chaotic dynamics have been observed. 

It is posible to use a simple theoretical scheme in order to explain qualitatively the transition of the system from 
non stationary regimes to stationary ones when fc is increased or Ly decreased. Such a scheme is based on a concept 
originated within the framework of nuclear reactor physics that simplifies noticeably the analysis of the so called 



4 



criticality problems or neutron density distributions, in many situations of interest. This is the concept of geometrical 
buckling. It helped to reduce the dimensionality of the problem, by assuming that directions transverse to the more 
relevant one can be essentially described by the fundamental mode. Such an assumption allows to include the effect 
of those extra dimensions through (lateral) "sink" contributions (that is taking into account the lateral leak) in a 
(typically one-dimensional) diffusion equation for the neutron density on the relevant direction IT^ . For instance, 
using this approximation it is possible to estimate the critical mass or radius of a nuclear system. Clearly, such an 
approximation gives reliable results due to the linearity of the neutron diffusion equation (at least over time scales 
shorter than those at which the transport coefficients -cross section, diffusion constants- change). 

Our system of equations is clearly nonlinear, hence such an approach could in principle not be applied. However, 
we qualitatively introduce such an idea forcing our equations to reduce their dimensionality and including associated 
sink or leak terms. 

For the system in Eq.(|^) in the rectangular geometry with Ly << L^, and albedo conditions at y = iLy given by 

dyu{x±, Ly) ~ ^ku{x, Ly), Dydyv{x,±Ly) — ^kv{x,Ly), (7) 

we will assume that the solution can be written as 

u{x, y, t) = u{x, t) cos((j„y) 

y, t) = v{x, t) cos{q^y). (8) 

Hence, we are approximating the profiles of the distributions of u and v in the y coordinate as being independent of 
X and t except for a global factor. The cosine functions give a more or less desirable form for the densities, having 
a maxima at y = 0. Clearly, the {k = 0)-case corresponds to g„ = q„ = with solutions independent of y. We may 
expect that our approximation will give qualitatively better results for small k. From Eq. (|^) we find the following 
relation that determines g„ and as function of k, Ly and D^: 

k ^ Qu tan{quLy) k = qyDy tan{qyLy), (9) 

and from Eq.(H) we obtain 

u 

V 

Now, since Ly << L^, 'we might think that, whether the dynamics is stationary or not will depend only on k and Ly 
and not on what happens on the x coordinate JTst , and we may consider the system as independent of x. We then 
obtain the following equations for a space independent activator-inhibitor system: 

" = /(") - qlu - V 

V ^ u — Dyq^v — "fv — c. (11) 

This system, for the chosen values of /(u), 7, and c, is oscillatory for null or small values of qu and qy, but there exist 
critical values of these quantities (which depend on Ly, Dy and k through Eqs.(^) beyond which the system begins 
to be monostable, i.e. asymptotically stationary. We associate this stationarity to that of our extended system. In 
Fig. 8 we have plotted with a dotted line the limit from stationary to non stationary regimes predicted in this way. 
Although the results do not agree quantitatively with those coming from simulations in the extended bidimensional 
system (solid line) , the occurrence of the phenomenon is clearly well predicted qualitativelly. 

IV. FINAL REMARKS 

We have analyzed the influence of partially absorptive boundary conditions in a typical reaction diffusion-model of 
an oscillatory active medium with different geometries. We have observed that the homogeneous periodic behaviour 
expected for the case of Neumann boundary conditions, when homogeneous initial condition are considered, is com- 
pletely altered when the absorption parameter is increased from zero. For small absorption, the behaviour continues 
being periodic although the homogeneity is lost and traveling waves are observed. For higher absorption (with pre- 
cise values of fc's depending on the diffusion constants), the periodicity is broken and quasiperiodic inhomogeneous 
oscillation and spatiotemporal chaos arise. 



+ f{u) 

- DyqyV + u 



"fV — c. 



(10) 



5 



The kind of phenomena described here is expected to appear in real systems sharing the properties of the model. 
The predictions made for the three systems analyzed (one-dimensional, circular and band-shaped) illustrates what 
can be expected in more complicated situations. For example, let us consider the same dynamics given by Eqs. 

with homogeneous initial conditions, but now in a square reactor with Neumann bc's everywhere except in a 
sector where albedo bc's are imposed. In Fig. 9 we show the contour plots of the u-fields for this system at several 
times taken in the asymptotic regime, assuming standard values of the parameters Dy and k. The spatiotemporal 
complexity arising is apparent. 

When analyzing the transition lines from periodic to nonperiodic regimes for the cases {k^ = k,ky — 0) and 
(ky = 0, fct, = k) we have observed shifts of these lines relative to the case fc„ = ky — k that are similar in both the 
one dimensional and the circular cases. These results imply that the system is more sensitive to boundary conditions 
on the w-field than on the u-field. 

As mentioned in the introduction, in a real experiment on oscillatory media, the waves emerging from the boundaries 
with imperfect reflectivity will interact with other patterns that may be spontaneously generated from fluctuations in 
the distribution of reactants in any part of the system. In this paper we have analyzed simpler situations that allowed 
us to isolate clearly the effects of the boundaries. The results obtained constitute a starting point for the study of 
more realistic situations. 
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FIG. 1. Propagation of a freezing front leading to the formation of a Turing pattern. Spatiotemporal plot of the u-field for 
L = 70, D„ = 2.4 and k„ = = 1. 

FIG. 2. Phase diagram ((_D„, fc)-plane) indicating the asymptotic behaviour of the one dimensional system for L = 70 and 
regions corresponding to inhomogeneous periodic oscillations, Turing patterns and quasiperiodic and chaotic 
behaviour, are indicated with PO, TP, QP, CH respectively. The transitions from periodic to nonperiodic regimes in the cases 
(fc„ = fc, fc„ = 0) and {ku = 0,kv = k) are indicated respectively with dashed and dotted lines. 

FIG. 3. Spatiotemporal plots of the u-fields in the asymptotic regime of the PO region. Calculations for 
D„ = 1.8, ku = fe„ = .02. 
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FIG. 4. Typical plots oi p{t) in for: (a) PO region = 1.8, ku = kv = .02); (b) QP region {D^ = 1., ku = ky = 1.); (c) CH 
region {Dy = 1.8, ku = kv = 1.). 

FIG. 5. Dynamics in the CH region. Top: spatiotemporal plot of the w-field for = 1.8, A;^ = /c„ = 1 in a temporal window 
of the asymptotic regime. Bottom: plot of p{t) for the same system in the same time window. 

FIG. 6. Phase diagrams indicating the asymptotic regimes of circular systems for: (a) L = 70; (b) L = 80. In both cases, 
the regions are indicated with the same nomenclature as in figure 2. With dashed and dotted lines we indicate the transitions 
from periodic to nonperiodic regimes in the cases {ku = k,kv = 0) and {ku = 0,kv = k) respectively. 



FIG. 7. Sketch of the band shaped system indicating the different boundary condition considered. 

FIG. 8. Diagram indicating the different behaviours observed in the band shaped system as function of k and Ly. The 
solid line indicates the transition between stationary and nonstationary regimes calculated from simulations. The dotted line 

corresponds to the estimation of this frontier with the discussed theoretical arguments. Crosses, squares, circles and triangles 
correspond to the observations in simulations of stationary, periodic, quasiperiodic and chaotic behaviour respectively. 

FIG. 9. Contour plots of the w-fields for a square system with Neumann bc's except in the right bottom corner (indicated 

with arrows in the left figure) where we have fixed albedo bc's with ku = kv = 1. From left to right the plots correspond to 
times t = 2000, 2004, 2008, 2012 and 2016. The other parameters are D„ = 2, and = Ly = 60. 
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FIG 1. Bouzat and Wio. 
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FIG 3. Bouzat and Wio. 





FIG 4. Bouzat and Wio. 
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